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Abstract 

We study numerically the spectrum and eigenfunctions of the quan- 
tum Neumann model, illustrating some general properties of a non trivial 
integrable model. 
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1 Introduction. 



The Neumann model, one of the first classical models to be shown integrable 
in a non trivial manner [Q, consists of a point constrained to live on a sphere, 
and submitted to different harmonic forces along the coordinate axis. It has 
been further studied, generalized to higher dimension and related to the Jacobi 
theory of geodesic flows on the ellipsoid by J. Moser 2^, while K. Uhlenbeck [3] 
found the conserved quantities in involution as required by Liouville theory. 
The classical Neumann model has been the object of much interest since its 
spectral curve is the generic hyperelliptic Riemann surface. Mumford ^ has 
been able to provide nice derivations of classical results on hyperelliptic curves 
by elaborating on this model and to give explicit formulae for its solution. 

Quantization of this model is non-trivial due to the constraints on the co- 
ordinates. This problem has been solved and the commutation of the natural 
quantization of the conserved quantities has been proved in This result 
was not obvious, but can be understood naturally at present in the framework 
presented in [H], which asserts that under fairly general conditions, classical 
integrable systems have integrable quantization. 

The common spectrum of the conserved quantities has been studied in [Jj, 
where it was shown that the Schrodinger equation separates using Neumann's 
coordinates. The problem reduces to the solution of a single differential equation 
involving the common set of eigenvalues, and one has to impose univaluedness of 
the wave function on the sphere to fix the eigenvalues. This is a difficult problem 
since, although the differential equation looks like a simple generalization of the 
hypergeometric equation, it cannot be solved by generalizations of Riemann's 
integral formulae, hence one cannot compute the monodromies that would be 
necessary to state the eigenvalue equations in closed form. However, writing 
the semi-classical Bohr-Sommerfeld conditions is manageable, and has been 
achieved in [7j. This line of approach has been developed further by D. Gurarie 
in [Hj. Introducing the Maslov indices for the classical trajectories, he has been 
able to refine the semi-classical quantization condition. He further considers 
the perturbative expansion around the trivial case where all oscillator forces 
vanish, i.e., one deals with free motion on a sphere. 

However, all these works fell short of giving concrete solutions for the spec- 
trum of the quantum Neumann model. In fact, almost no non-trivial quantum 
integrable model yields exact formulae for the spectrum, the Calogero model 
being a notable exception. The aim of this paper is to present a fairly detailed 
study of the spectrum and eigenstates of the quantum Neumann model, ob- 
tained through numerical methods. The problem has been made manageable 
by the recognition that, for the model on S^~^, N parity operators commute 
with all conserved quantities. This translates in a unique set of boundary con- 
ditions for the factors in the separated wave function, so that these factors are 
slices of a unique solution of the separated Schrodinger equation. A local study 
allows to translate the numerical problem to a multiparameter spectral prob- 
lem with an assorted number of boundary conditions, for which the COLNEW 
program yields a very efficient solution. 

We shall first recall the definition of the classical Neumann model, introduc- 
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ing a two by two Lax pair formulation. This allows to introduce in a natural 
way Neumann's separating coordinates and Uhlenbeck's conserved quantities. 
We then describe the formulation of the quantum problem, followed by a pre- 
sentation of the numerical methods. Finally, we show the first thirty-six energy 
levels of the system plotted as a function of the strength of the potential. 



2 The Neumann model. 



We begin by presenting the classical Neumann model, which is discussed in a 
different manner in the book 0. We start from a 2A^-dimensional free phase 
space {xn, Un, n = 1- ■ ■ N} with canonical Poisson brackets: {xn,ym} = ^nm 
and introduce the "angular momentum" antisymmetric matrix: J^i = x^Ui — 
xiUk and the Hamiltonian: 

^ = \Y.'^ki + lY.^'kxi (1) 

kytl k 

We shall assume in the following that: ai < a2 < ■ ■ ■ < a^- With the vectors 
X = (xk) and Y = (yi^) and the diagonal constant matrix A = (a^dki), the 
Hamiltonian equations read: 

X = -JX Y = -JY-AX 

They automatically ensure that = Y1 

x| remains constant and lead to the 
non-linear Newton equations for the particle: 

Xk = -r^akXk - Xk ^(x/^/r^ - aixf) 
I 

The Liouville integrability of this system is a consequence of the existence 
of (A^ — 1) independent quantities in involution 0: 

F. = xi + E^ (2) 

This can be understood easily when one introduces the following Lax pair for- 
mulation of the Neumann model: 

iW = i» + EF^i- io=(_"i °), L* = »(2 (3) 

As usual in such a formulation, all dynamical variables are described by the gk- 
Since L^. is nilpotent, can be taken as living in the group of lower triangular 
matrices of unit determinant and we parametrize it as: 

9^=(? (4) 



Vk 

The matrix L{t) can therefore be written: 
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with 



k ^ ~ "'^ k ^ ~ '^^ k ^ ~ '^'^ 



The spectral curve, det(L(t) — si) = 0, reads 



„2 



v{tf + u{t)w{t)=0 = s^ + Yj^^ (6) 

k ^ 



This is the equation of a smooth hyperelhptic Riemann surface of genus g = N— 
1 We see that the Fk are integrals of motion since they appear in the moduli 
of the spectral curve. Notice that Fk = J2k and H = 1/2 o^Ffc. 

Since ^^i. x| commutes with the Hamiltonians, it generates a canonical 
transformation : 

{xk,yk) — > + 2Axfc) (7) 

which is a symmetry of the system. The Hamiltonian reduction with respect to 
this symmetry is the constrained system on the sphere. In the Lax setup, this 
symmetry corresponds to conjugating L{t) by the constant matrix 

1 
2A 1 

Note that this conjugation fixes Lq. The properly reduced system has a phase 
space of dimension 2(A^ — 1) = 2g. 

According to the standard Sklyanin procedure JOl , we introduce separating 
coordinates as the roots of Li2{t) = u{t). The positivity of the xf ensures that 
there is exactly one root in each of the intervals ]ak, afc+i[ so we can set: 

oi ^ ^1 ^ ^2 ^ ■ ■ • ^ o-N-i ^ tN-i < qn 

The coordinates tj in this domain form an orthogonal system which covers 
bijectively the quadrant Xi > 0, as illustrated in Fig. 1: 

V nfc(*-afc) V Ui^ki'^k - ai) 



When t ^ Uk we have Xk — — ak so that other quadrants are explored by 
analytic continuation around the branch point at t = ak- 

3 The Schrodinger equation. 

We shall seek wave functions which factorize in the separating variables tk- 

^{h, • • • , tiv-l) = ^l(il) ^2(^2) • • • ^N^l{tN-l) (9) 

the Hilbert space being the space of functions of variables ti with the measure: 



■Jldti 
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Figure 1: A uniform grid in the separating coordinates for = 3. 

It has been shown in [Jj that in order to simultaneously solve the eigenvalue 
equations i^fc^ = fk^ all the satisfy the same one-dimensional Schrodinger 
equation: 



This is a linear differential equation with singularities at the Uk and at infinity. 
More precisely at each there is a regular singularity jllj . with exponents 
and 1/2, while the singularity at infinity is irregular. Hence this equation 
appears as a generalization of the hyper geometric equation, however unlike the 
case of three regular singularities there doesn't appear to exist integral formulae 
for the solution, precluding to perform explicit analytic continuation around the 
branch points at ak- 

There exists, fortunately, a workaround. Remark that all Hamiltonians Fk 
are invariant under the N parity transformations {xj,yj) — > {—Xj,—yj) so that 
one can simultaneously diagonalize both the Hamiltonians and the parities. 
That is, we can require that the wave function is even or odd in Xj for each j, or 
equivalently that is is a function of the x| multiplied by a monomial Y[ with 
rij S {0, 1}. Let us now recall that if to is a regular singularity with exponents 
a, f3, there exists convergent developments of the solution of the form 



where m is a or /?. The exponents being 0, 1/2 at each a^, we see that these 
solutions are respectively even or odd under —Xk, since x^ ~ y/t — for 

t Ofc. 




(10) 



{t - toTipQ +pi{t - to) +P2{t - tof + ■■■) 
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Moreover let us look at the parity properties of the wave function ^ in eq.® 
with respect to Xk- There are two ways to approach Xk = 0: 

tk — *■ «fc or Ofc Qk-i < tfe-i < <tk < ct/fc+i 

Requiring that the product ^ k-i{tk-i)^ k{tk) leads to a definite parity state in 
Xk amounts to saying that ^k-iit) ^-nd ^'^(t) share the same exponent at a^, 
hence both share the same analytic expansion in the interval ]afc-i,afc+i[- By 
induction, we arrive at the stronger statement that in order to have definite 
parities with respect to the x^, all functions ^'fc in eq.® are in fact one and 
the same solution ^{t) of the differential equation on the whole interval 
]ai,aAr[, but constrained to have definite exponents at the points at-, so as to 
produce the above monomial W x"^ . This is the constraint which quantizes 
the -Ffc in ea. (|lUI) . First it is clear that a wave function ^ so obtained on the 
quadrant X/t > extends by parity to a univalued wave function on the sphere 
satisfying the Schrodinger equation. Second, start from a solution with definite 
exponent at ai. In general its continuation to 02 will be a superposition of the 
two pure exponent solutions. Requiring that there is no mixing imposes one 
condition on the fk- Similarly at 03, . . . , cat so that we get A'" — 1 conditions on 
the — 1 independent f^, the quantization conditions. 

It remains to express these conditions in a form suitable for numerical anal- 
ysis. We can always reduce the problem to the case where the solution is even 
under all Xk —Xk- Indeed if we want an odd solution, we have only to write 
^{t) = \Jt — ak^{t) and obtain the differential equation satisfied by which 
is very similar to eg. (110(1 . Then we require that $ be even. Hence we need to 
express that the solution "^{t) has a regular power series development at each 
point afc. Inserting ^'(t) = Po + pi{t — to) + p2{t — to)^ + • • • in ea. ((TU)) and 
requiring that the polar terms cancel we get the equation pi = pQfk/{2h^), that 
is: 

^'{at) = ^{ak)fk/{2h^) (11) 

Of course these equations and ea. (|lU() are linear in ^' so a definite solution for 
\I' is obtained when requiring for example a normalization condition such as 
^(ai) = 1. 

4 Numerical techniques. 

The problem is now set in a form which happens to be directly tractable by 
a known numerical analysis package called COLNEW^. This package allows to 
solve so-called mixed order boundary values problems: suppose we have a set 
of differential equations of various orders u^"^^^ = fi{u) where u is a vector of 
all functions Ui and all their derivatives up to order — 1 for Uj, where / 
may be linear or non linear. The general solution depends on ^ rrij constants. 
The package allows to impose ^ mj "boundary conditions" and determines the 
solution satisfying all of them. These boundary conditions have to be of the 

^This is work of Ascher, Christiansen, Russell, and Bader, which can be found at 
|http: //www. netlib.org/ode/, 
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form gj[u){C,j) = where Qj are expressions of the same vector u as above, and 
the points Cj are ^ nrii points arbitrarily disposed in the considered domain 
or its boundary. The trick to apply this scheme to our problem is to treat 
the quantities fk as functions of t which happen to be constant, i.e. satisfy 

= 0. Supplementing ea. (|lUjl by these — 1 independent conditions we 
get a mixed order non linear differential system with total order 2 + {N — 
l).l = + 1. The boundary conditions we impose are the N above equations 
^'{cLk) = ^[o,k)fk/{'^h'^) at each plus a normalization such as ^(ai) = 1 
which completely fixes things. 

Finally to appreciate how the numerical package can deal with the singu- 
larities at the Qk appearing in eq. (|inj) . it is in order to say a few words about 
it. The computation starts by cutting the considered interval into a number of 
subintervals. One is free to require that the are included in the subdivision. 
In each interval, the solution is approximated by some polynomial. The coeffi- 
cients of these polynomials are fixed by a set of equations: first one requires that 
the differential equation is exactly satisfied at a number of so-called "colloca- 
tion points" inside the intervals. These points are chosen at Gaussian positions 
as in numerical integration, in order to improve convergence. Note that, these 
points being interior while the are boundary points, we are always away from 
the poles in l/(t — Ofc) when computing these equations. Second the package 
imposes that the solutions from interval to interval fit in a smooth way, and 
that the boundary conditions are obeyed. These conditions don't suffer from 
divergent factors at and moreover imply in a finite way that the differential 
equation is indeed obeyed at the a/, and that the solution has the appropriate 
parity property. 

Things are set up so that one gets a non linear system for the coefficients, 
with as many equations as unknowns. This system is then solved by Newton- 
Raphson method, an error is evaluated, and a more refined solution is computed 
until a defined tolerance is obtained on all the equations. In practice we have 
found that the package works remarkably well for our problem, and is able to 
compute a fair number of solutions in a couple of seconds on a modern machine. 
Moreover its use is enhanced by the existence of a wrapper in Scilab^ which 
allows to program these computations in a convenient way. 

The numerical results below will be restricted to the case A'^ = 3. There 
would not be any difficulty to study similarly higher dimensional problems. To 
reduce the parameter space to explore, we first normalize coefficients. Translat- 
ing all the Ofc by the same quantity does not change the and adds a constant 
to the energy. We therefore set ai = 0. The equation is also invariant by 
a common rescaling of t and the by A and the rescaling of the fk by 1/A. 
The h factor can also be absorbed in the fk so that we are reduced to two 
parameters. We choose 02 = 1 and have a parameter y = 03/02 characterizing 
the asymmetry of the model. Now Ylk fk ~ '^^ where 




(12) 



k 



^http: / /scilabsoft .inria.fr / 
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This single parameter can either be seen as characterizing the strength of the 
potential energy compared to the kinetic energy, or the approach to the semi- 
classical regime, which are therefore equivalent in this case. Finally, let us write 
the equations to be solved: 

(P 1/1 1 1 \ d 

+ ^ - + -. 7+ ^ 



dt"^ 2\t t-1 t-yJdt 



t t-1 t-y 



'3 1 I \ d 

+ ^\- + -. 7 + 



dt'^ 2\t t-1 t-yJdt 

-f2-h + l/^ + l/m , /2^-V4 , /3-V(%) )yp(,)_o (14) 

The second equation is the one we obtain when setting ^ ^/t"^. There 
exist six similar equations describing, with the two displayed, the eight possible 
combinations of parities of the solution. Equation (|13|) is known as Wangerin's 
equation, but its properties have not been explored. When v = the singularity 
at infinity becomes regular, and the equation reduces to Lame's equation. A 
great deal of work on this equation is summarized in and is relevant for 
our study as we see below. 



5 Numerical results. 

A single spectrum is not really instructive. We shall therefore show the evolution 
of the 36 lowest energy levels when the parameter v grows. The number of 
different energy levels displayed is limited for the readability of the plots. States 
are characterized by their sets of parities and by the number of zeros of the wave 
function in each of the two intervals ]0, 1[ and ]l,y[. The principal difficulty is 
to be able to choose a particular solution with a given number of zeros. The 
COLNEW package allows for a powerful mean to achieve this: an initial guess of 
the solution can be provided. However good guesses are not easy to find except 
for V = 0. In this case, the theory of the Lame spheroidal harmonics teaches 
us that the solutions which are analytic at the ak are polynomials. The zeros 
of these polynomials satisfy simple equations and it is therefore not difficult to 
give numerical approximations of the solutions. To reach higher values of w, we 
have used the continuation method, which consists in taking the solution for 
the old value of v as the initial guess for the next value of v. If the difference 
between successive values of v is sufficiently small, the new solution has the 
same quantum numbers^. Below we have plotted E/2 = /2 + y/3 but we could 
have plotted individually each function of v. 

To have an idea of the effect of the asymmetry parameter y, we plot the 
energies for three values of this parameter. The line styles correspond to the 
parity properties of the states. The heavy solid lines correspond to the states 

^The standard COLNEW wrapper in Scilab is however unsuitable to perform this contin- 
uation and we had to modify it. 
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Figure 2: Energy as a function of v. Case y = 1.1. 



which are completely even, the other hnes are sohd or dashed according to 
the total parity of the states. In Fig. 2, y is equal to 1.1, meaning that the 
two coordinates X2 and are nearly equivalent. In fact, when = a^+i 
there is a singular term in and F^+i, but we can change for their sum 
and difference, with the sum being regular and the difference proportional to 
Jkk+v ^^^^ case of enhanced symmetry, the algebro-geometric solution of 
the classical model becomes singular, see jH], but the quantum solution has no 
particular trouble. Note however that a degeneracy of energy levels is apparent 
on the figure. The case y = 2 in Fig. 3 is representative of the case where 
the oscillator strengths are equally distributed, while the case y = 5 in Fig. 4 
illustrates the case where one of them is bigger. A remarkable feature of these 
plots is that the levels converge at v = to integer or half integer values. 
This is easily understood: = is equivalent to setting the potential to zero, 
hence the problem has spherical symmetry, so the energy is equal to a factor 
times the eigenvalue of the operator, j{j + 1) with a 2j + 1 degeneracy. 
This is easily observed on the curves. An other feature of these plots is that 
pair of levels converge for large v. This can be understood from the potential 
barrier at the xi = plane, so that the probability is concentrated around 
the opposite poles. In this limit, there are therefore two equivalent states, just 
separated by the exponentially small tunneling probability. When v becomes 
even greater, the states occupy but a small patch around the poles and the effect 
of the curvature of the sphere becomes small. The problem reduces to the one 
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Figure 3: Energy as a function of v. Case y = 2. 

of two independent harmonic oscillators, which have energies proportional to 
^yv. Remark however that with growing quantum numbers, the level crossings 
between the two regimes v small and v large occur at larger values of f , so that 
for a given value of v, the eigenfunctions of sufficient energies are perturbations 
of the spheroidal harmonics. 

Finally the solution of eq. (|lfl|) for an excited state with v = 2 and y = 2 
is plotted in figure ©• Precisely for this solution, which has one zero in ]0, 1[ 
and two zeroes in ]l,y[ we have /2 = —0.338, /s = 9.913 hence the "energy" 
is /2 + yfs = 19.49. Recall that the wave function is ^{ti)^{t2), with < 
ti < 1 < t2 < y. Note that the locus of zeroes of the wave function is formed 
of intersecting lines due to the separation of variables. This is clearly non 
generic and would be destroyed by any perturbation. The existence of such 
singular node curves in the eigenstates of a quantum Hamiltonian is evidence 
for integrability. 

6 Conclusion. 

An efficient calculation of the spectrum of the quantum integrable Neumann 
model has been displayed. A number of degeneracies are evident from the 
above figures. It will be the purpose of a separate work to explicit these phe- 
nomena and relate them with the semiclassical analysis of the equation. This 
makes contact with the solution of the classical case in terms of hyperelliptic 
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functions. The interesting features are the transitions between different semi- 
classical quantization conditions. 
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